{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "82b847e6d9e2cc0",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:26.978496Z",
     "start_time": "2024-07-24T08:59:26.419317Z"
    }
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import blosc2 as b2\n",
    "import numpy as np\n",
    "import tables as tb\n",
    "import os\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a80c910ca02f49e0",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:26.982122Z",
     "start_time": "2024-07-24T08:59:26.979712Z"
    }
   },
   "outputs": [],
   "source": [
    "# A tomography-like array: a stack of 2D images (greyscale).\n",
    "# Each image corresponds to a chunk in the array.\n",
    "shape = (10, 25600, 19200)\n",
    "dtype = np.dtype(\"u2\")\n",
    "# shape = (100, 256, 256)  # for tests\n",
    "chunkshape = (1, *shape[1:])\n",
    "clevel = 5\n",
    "# The next can be set as environment variables or directly here.\n",
    "# To select the compressor, and compression level for the regular HDF5 Blosc2 filter\n",
    "os.environ[\"BLOSC_CLEVEL\"] = str(clevel)\n",
    "os.environ[\"BLOSC_COMPRESSOR\"] = \"zstd\"\n",
    "# You can play with different filters here (NOSHUFFLE/SHUFFLE/BITSHUFFLE)\n",
    "# os.environ[\"BLOSC_SHUFFLE\"] = \"SHUFFLE\"\n",
    "# You can set the number of threads for Blosc\n",
    "# os.environ['BLOSC_NTHREADS'] = str(b2.nthreads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c7307a8bb88bddd8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:26.984387Z",
     "start_time": "2024-07-24T08:59:26.982762Z"
    }
   },
   "outputs": [],
   "source": [
    "# Blosc2 blockshape cannot be specified via `tb.Filters` because the interface is\n",
    "# not ready for passing it yet, so to compare apples with apples, avoid using it.\n",
    "# However, you may want to experiment with different blockshapes with direct chunking.\n",
    "# b2_blockshape = (1, *tuple(d // 2 for d in chunkshape[1:]))  # 4 blocks per chunk\n",
    "# b2_blockshape = (1, 800, 600)  # around 1 MB per block\n",
    "# blocksize = np.prod(b2_blockshape) * dtype.itemsize\n",
    "# print(\"Blocksize:\", blocksize, b2_blockshape)\n",
    "# Use the next line to set the blocksize for regular I/O.  Unfortunately, this is not\n",
    "# compatible with a multidim blockshape in direct chunking (we are sending different\n",
    "# requirements to the Blosc2 library), so better don't use it.\n",
    "# os.environ['BLOSC_BLOCKSIZE'] = str(blocksize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7e7d2e7bcfcd464a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:27.321321Z",
     "start_time": "2024-07-24T08:59:26.985740Z"
    }
   },
   "outputs": [],
   "source": [
    "# The 'image' dataset that will be (repeatedly) stored in the HDF5 file.\n",
    "np_data = np.arange(np.prod(chunkshape), dtype=dtype).reshape(chunkshape)\n",
    "\n",
    "# The measured speeds\n",
    "speeds = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2073fbed989f97bc",
   "metadata": {},
   "source": [
    "### Direct chunking\n",
    "\n",
    "Let's start by creating a new HDF5 file using direct chunking.  The procedure is as follows:\n",
    "1. Create an extendable array.\n",
    "2. Grow the array to the final size.\n",
    "3. Create the compressed data chunk (cframe).\n",
    "4. Write the compressed data chunk by chunk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "337bbc38005ef881",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:32.757711Z",
     "start_time": "2024-07-24T08:59:27.341473Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[direct chunking] Wrote 10 chunks in (5.71 s). Speed: 1722.51 MB/s.\n"
     ]
    }
   ],
   "source": [
    "fname = \"direct-chunking.h5\"\n",
    "with tb.open_file(fname, mode=\"w\") as h5f:\n",
    "    array = h5f.create_earray(\n",
    "        \"/\",\n",
    "        \"array\",\n",
    "        atom=tb.Atom.from_dtype(dtype),\n",
    "        shape=(0, *shape[1:]),\n",
    "        # Setting both args tells others that data is compressed using Blosc2,\n",
    "        # and it should not be handled as plain data.\n",
    "        filters=tb.Filters(complevel=clevel, complib=\"blosc2\"),\n",
    "        chunkshape=chunkshape,\n",
    "    )\n",
    "\n",
    "    # First, grow the array without actually storing data.\n",
    "    array.truncate(shape[0])\n",
    "    # Now, do store the data as raw chunks.\n",
    "    coords_tail = (0,) * (len(shape) - 1)\n",
    "\n",
    "    def do_write():\n",
    "        for c in range(shape[0]):\n",
    "            # The same image/chunk will be written over and over again.\n",
    "            b2_data = b2.asarray(\n",
    "                np_data, chunks=chunkshape, cparams=dict(clevel=clevel)\n",
    "            )\n",
    "            wchunk = b2_data.to_cframe()\n",
    "            chunk_coords = (c,) + coords_tail\n",
    "            array.write_chunk(chunk_coords, wchunk)\n",
    "\n",
    "    start = time.time()\n",
    "    do_write()  # cProfile.run('do_write()')\n",
    "    elapsed = time.time() - start\n",
    "    speed = np.prod(shape) * dtype.itemsize / elapsed / 1e6\n",
    "    print(\n",
    "        f\"[direct chunking] Wrote {shape[0]} chunks in ({elapsed:.3} s).\"\n",
    "        f\" Speed: {speed:.2f} MB/s.\"\n",
    "    )\n",
    "    speeds[\"direct\"] = {}\n",
    "    speeds[\"direct\"][\"write\"] = speed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e123e547a5cd3d1a",
   "metadata": {},
   "source": [
    "Now, let's read the data back.\n",
    "\n",
    "The procedure is as follows:\n",
    "\n",
    "1. Read the compressed data chunk.\n",
    "2. Decompress the data chunk.\n",
    "3. Read the data into a numpy array from the decompressed data chunk.\n",
    "4. Repeat for all chunks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ff10733c19ef2c8f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:35.666059Z",
     "start_time": "2024-07-24T08:59:32.759037Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[direct chunking] Read 10 chunks of 9375.0 MB (1.9 s). Speed: 5172.77 MB/s.\n"
     ]
    }
   ],
   "source": [
    "with tb.open_file(fname, mode=\"r\") as h5f:\n",
    "    array = h5f.root.array\n",
    "\n",
    "    coords_tail = (0,) * (len(shape) - 1)\n",
    "\n",
    "    def do_read():\n",
    "        tsize = 0\n",
    "        for c in range(shape[0]):\n",
    "            chunk_coords = (c,) + coords_tail\n",
    "            rchunk = array.read_chunk(chunk_coords)\n",
    "            ndarr = b2.ndarray_from_cframe(rchunk)\n",
    "            np_data2 = ndarr[:]\n",
    "            tsize += np_data2.size\n",
    "        return tsize * dtype.itemsize\n",
    "\n",
    "    start = time.time()\n",
    "    tsize = do_read()\n",
    "    elapsed = time.time() - start\n",
    "    assert tsize == np.prod(shape) * dtype.itemsize\n",
    "    speed = np.prod(shape) * dtype.itemsize / elapsed / 1e6\n",
    "    print(\n",
    "        f\"[direct chunking] Read {shape[0]} chunks of {tsize/2**20} MB ({elapsed:.3} s).\"\n",
    "        f\" Speed: {speed:.2f} MB/s.\"\n",
    "    )\n",
    "    speeds[\"direct\"][\"read\"] = speed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9a04e0d34d85659",
   "metadata": {},
   "source": [
    "## Regular chunking\n",
    "\n",
    "Let's now create a new HDF5 file using regular chunking.  The procedure is as follows:\n",
    "1. Create an extendable array, using Blosc2 as a compressor.\n",
    "2. Grow the array to the final size.\n",
    "3. Write the data chunk by chunk (the HDF5 pipeline will call Blosc2 automatically)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f60cdacf9a59b39c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:51.756959Z",
     "start_time": "2024-07-24T08:59:35.667420Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[regular chunking] Wrote 10 chunks in (21.5 s). Speed: 458.24 MB/s.\n"
     ]
    }
   ],
   "source": [
    "fname2 = \"regular-chunking.h5\"\n",
    "with tb.open_file(fname2, mode=\"w\") as h5f:\n",
    "    array = h5f.create_earray(\n",
    "        \"/\",\n",
    "        \"array\",\n",
    "        atom=tb.Atom.from_dtype(dtype),\n",
    "        shape=(0, *shape[1:]),\n",
    "        # Setting both args tells others that data is compressed using Blosc2,\n",
    "        # and it should not be handled as plain data.\n",
    "        filters=tb.Filters(complevel=clevel, complib=\"blosc2\"),\n",
    "        chunkshape=chunkshape,\n",
    "    )\n",
    "\n",
    "    # First, grow the array without actually storing data.\n",
    "    array.truncate(shape[0])\n",
    "    # Now, do store the data as raw chunks.\n",
    "    coords_tail = (0,) * (len(shape) - 1)\n",
    "\n",
    "    def do_write():\n",
    "        for c in range(shape[0]):\n",
    "            # The same image/chunk will be written over and over again.\n",
    "            array[c] = np_data\n",
    "\n",
    "    start = time.time()\n",
    "    do_write()  # cProfile.run('do_write()')\n",
    "    elapsed = time.time() - start\n",
    "    speed = np.prod(shape) * dtype.itemsize / elapsed / 1e6\n",
    "    print(\n",
    "        f\"[regular chunking] Wrote {shape[0]} chunks in ({elapsed:.3} s).\"\n",
    "        f\" Speed: {speed:.2f} MB/s.\"\n",
    "    )\n",
    "    speeds[\"regular\"] = {}\n",
    "    speeds[\"regular\"][\"write\"] = speed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b7cf7c9cc2fca0b",
   "metadata": {},
   "source": [
    "Now, let's read the data back.\n",
    "\n",
    "The procedure is as follows:\n",
    "\n",
    "1. Read the compressed data chunk. The HDF5 pipeline will call Blosc2 automatically for decompression.\n",
    "2. Repeat for all chunks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "522b828c5c889de9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:55.397480Z",
     "start_time": "2024-07-24T08:59:51.762261Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[regular chunking] Read 10 chunks of 9375.0 MB (8.35 s). Speed: 1176.62 MB/s.\n"
     ]
    }
   ],
   "source": [
    "with tb.open_file(fname2, mode=\"r\") as h5f:\n",
    "    array = h5f.root.array\n",
    "\n",
    "    coords_tail = (0,) * (len(shape) - 1)\n",
    "\n",
    "    def do_read():\n",
    "        tsize = 0\n",
    "        for c in range(shape[0]):\n",
    "            np_data2 = array[c]\n",
    "            tsize += np_data2.size\n",
    "        return tsize * dtype.itemsize\n",
    "\n",
    "    start = time.time()\n",
    "    tsize = do_read()  # cProfile.run('do_read()')\n",
    "    elapsed = time.time() - start\n",
    "    assert tsize == np.prod(shape) * dtype.itemsize\n",
    "    speed = np.prod(shape) * dtype.itemsize / elapsed / 1e6\n",
    "    print(\n",
    "        f\"[regular chunking] Read {shape[0]} chunks of {tsize/2**20} MB ({elapsed:.3} s).\"\n",
    "        f\" Speed: {tsize / elapsed / 1e6:.2f} MB/s.\"\n",
    "    )\n",
    "    speeds[\"regular\"][\"read\"] = speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2b0be52ed78882c0",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-24T08:59:55.541373Z",
     "start_time": "2024-07-24T08:59:55.399537Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGzCAYAAADOnwhmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEVElEQVR4nO3dd3wUdf7H8feGJEtIJSQklAABlCJFygG5E4GjRAQUxBMbTRCBoAZORDyl6CnKWcCCnnIS60k7FVDAUAWJikAQkC4dEhBIoSWQfH9/8Mj8soaShYQNmdfz8ZjHg5357sxndnYyb2a+M+swxhgBAADYmJenCwAAAPA0AhEAALA9AhEAALA9AhEAALA9AhEAALA9AhEAALA9AhEAALA9AhEAALA9AhEAALA9AhHgQePGjZPD4fB0GSikZcuWyeFwaNmyZZ4uRZKUkJAgh8Ohn3/+2SPLr1Gjhrp27XrZdg6HQ+PGjSv+goCrQCACikjewSlvKFu2rCpXrqzY2Fi98cYbyszM9HSJllOnTmncuHEl7sCeN3h7e6tKlSrq16+fDhw44OnyANiAt6cLAEqb5557TtHR0Tp79qxSUlK0bNkyxcfH67XXXtOcOXPUqFEjq+0zzzyjp5566prXeOrUKY0fP16S1LZt22u+/IvJ++zOnDmjH374QQkJCVq5cqU2btyosmXLero8XKHTp0/L25vDDUo2vqFAEevcubOaN29uvR49erSWLFmirl276o477tDmzZvl5+cnSfL29r7sgSI3N1fZ2dm2CAT5P7uBAwcqLCxML7/8subMmaN77rnHw9UVvZMnT8rf39/TZRQ7O3x3cf3jkhlwDfz1r3/Vs88+qz179uiTTz6xxl+oD5HD4dCwYcP06aef6qabbpLT6dSCBQskSQcOHNBDDz2kiIgIOZ1O3XTTTfrggw8KLO/MmTMaN26cbrzxRpUtW1aVKlXSXXfdpZ07d2r37t0KDw+XJI0fP966THWxPh4///yzHA6HPvzwwwLTFi5cKIfDoXnz5kmSMjMzFR8frxo1asjpdKpixYrq2LGj1q5de0WfW+vWrSVJO3fudBm/ZcsW3X333QoNDVXZsmXVvHlzzZkzp8D7f/nlF7Vp00Z+fn6qWrWq/vnPf2ratGlyOBzavXu31e5i61+jRg3169fvkjWuWLFCf/vb31StWjU5nU5FRUVp+PDhOn36tEu7fv36KSAgQDt37tTtt9+uwMBAPfDAA5ec94EDBzRgwABVrlxZTqdT0dHRGjJkiLKzs13aZWVlacSIEQoPD5e/v7969OihI0eOuLQp7DrmXb78/vvvLzvPC/nwww/l7e2tkSNHXnTZed/7HTt2qF+/fgoJCVFwcLD69++vU6dOuczv9OnTeuyxxxQWFqbAwEDdcccdOnDgAP2SUOQ4QwRcI71799bTTz+tb7/9Vg8//PAl2y5ZskQzZszQsGHDFBYWpho1aig1NVWtWrWyAlN4eLjmz5+vAQMGKCMjQ/Hx8ZKknJwcde3aVYsXL9a9996rxx9/XJmZmUpMTNTGjRvVoUMHvfPOOxoyZIh69Oihu+66S5JcLuXl17x5c9WsWVMzZsxQ3759XaZNnz5d5cuXV2xsrCRp8ODBmjVrloYNG6b69evr6NGjWrlypTZv3qymTZu6/ZnlhZby5ctb4zZt2qS//OUvqlKlip566in5+/trxowZ6t69u2bPnq0ePXpIOh8m2rVrJ4fDodGjR8vf319Tp06V0+l0u45LmTlzpk6dOqUhQ4aoQoUK+umnn/Tmm29q//79mjlzpkvbc+fOKTY2VrfccoteeeUVlStX7qLzPXjwoFq0aKG0tDQNGjRIdevW1YEDBzRr1iydOnVKvr6+VttHH31U5cuX19ixY7V7925NmjRJw4YN0/Tp0694va5knu+9954GDx6sp59+Wv/85z8vu4x77rlH0dHRmjBhgtauXaupU6eqYsWKevnll602/fr104wZM9S7d2+1atVKy5cvV5cuXa54vYCLMgCKxLRp04wks3r16ou2CQ4ONk2aNLFejx071vxxN5RkvLy8zKZNm1zGDxgwwFSqVMn8/vvvLuPvvfdeExwcbE6dOmWMMeaDDz4wksxrr71WYPm5ubnGGGOOHDliJJmxY8cWat1Gjx5tfHx8zLFjx6xxWVlZJiQkxDz00EMu6xcXF1eoeeaX99ktWrTIHDlyxOzbt8/MmjXLhIeHG6fTafbt22e1bd++vWnYsKE5c+aMy3r9+c9/NjfccIM17tFHHzUOh8OsW7fOGnf06FETGhpqJJldu3ZZ4y/2WVSvXt307dvXer106VIjySxdutQal/e55zdhwgTjcDjMnj17rHF9+/Y1ksxTTz1VqM+kT58+xsvL64Lfp7ztmPe5dejQwRpnjDHDhw83ZcqUMWlpaW6vozvzrF69uunSpYsxxpjJkycbh8Nhnn/++QLL+OOy8773+b87xhjTo0cPU6FCBev1mjVrjCQTHx/v0q5fv35ufX+BwuCSGXANBQQEFOpuszZt2qh+/frWa2OMZs+erW7duskYo99//90aYmNjlZ6ebl2Wmj17tsLCwvToo48WmO+V3uLfq1cvnT17Vv/73/+scd9++63S0tLUq1cva1xISIh+/PFHHTx48IqW06FDB4WHhysqKkp33323/P39NWfOHFWtWlWSdOzYMS1ZskT33HOPMjMzrc/g6NGjio2N1fbt26270hYsWKCYmBjdfPPN1vxDQ0Mve5nKXXn9waTzfYJ+//13/fnPf5YxRuvWrSvQfsiQIZedZ25urr788kt169bNpT9anj9ux0GDBrmMa926tXJycrRnzx53VuWK5zlx4kQ9/vjjevnll/XMM88UehmDBw92ed26dWsdPXpUGRkZkmRdKh46dKhLuwt9t4GrRSACrqETJ04oMDDwsu2io6NdXh85ckRpaWl67733FB4e7jL0799fknT48GFJ5/vb1KlTp0jv6mncuLHq1q3rcrlk+vTpCgsL01//+ldr3MSJE7Vx40ZFRUWpRYsWGjdunH777bdCL+ftt99WYmKiZs2apdtvv12///67yyWuHTt2yBijZ599tsDnMHbsWEn//zns2bNHtWvXLrCMC427Gnv37lW/fv0UGhqqgIAAhYeHq02bNpKk9PR0l7be3t5WuLuUI0eOKCMjQw0aNChUDdWqVXN5nXeJ8fjx44V6/9XMc/ny5Ro1apRGjRrl0m+oKJaxZ88eeXl5FdgfinobAhJ9iIBrZv/+/UpPTy/UH/P8Zx2k82cMJOnBBx8s0I8nz8X6ABWVXr166YUXXtDvv/+uwMBAzZkzR/fdd59L8LrnnnvUunVrffHFF/r222/1r3/9Sy+//LL+97//qXPnzpddRosWLawzIt27d9ctt9yi+++/X1u3blVAQID1OTzxxBNWv6U/KsqDZU5OzmWnd+zYUceOHdOoUaNUt25d+fv768CBA+rXr59Vbx6n0ykvr6L/f2iZMmUuON4Yc9n3XmwdCzvPm266SWlpafr444/1yCOPFAgvl3I1dQNFjUAEXCMff/yxJF30QH4p4eHhCgwMVE5Ojjp06HDJtrVq1dKPP/6os2fPysfH54JtruTSWa9evTR+/HjNnj1bERERysjI0L333lugXaVKlTR06FANHTpUhw8fVtOmTfXCCy8UKhDlV6ZMGU2YMEHt2rXTW2+9paeeeko1a9aUJPn4+Fz2c6hevbp27NhRYPyFxpUvX15paWku47Kzs3Xo0KFLLmPDhg3atm2bPvzwQ/Xp08can5iYeMn3XU54eLiCgoK0cePGq5pPfle6jpcTFhamWbNm6ZZbblH79u21cuVKVa5c+armmad69erKzc3Vrl27dMMNN1jjL7QNgavFJTPgGliyZImef/55RUdHX1EfljJlyqhnz56aPXv2BQ+S+W+H7tmzp37//Xe99dZbBdrl/c877+6mPx4gL6VevXpq2LChpk+frunTp6tSpUq69dZbrek5OTkFLhFVrFhRlStXVlZWVqGXk1/btm3VokULTZo0SWfOnFHFihXVtm1b/fvf/77ggTz/5xAbG6ukpCQlJydb444dO6ZPP/20wPtq1aql7777zmXce++9d9kzRHlnOPKf0TDGaPLkyYVav4vx8vJS9+7dNXfu3Av+LMeVnEG50nUsjKpVq2rRokU6ffq0OnbsqKNHj171PKX//8/DlClTXMa/+eabRTJ/ID/OEAFFbP78+dqyZYvOnTun1NRULVmyRImJiapevbrmzJlzxQ+pe+mll7R06VK1bNlSDz/8sOrXr69jx45p7dq1WrRokY4dOyZJ6tOnjz766CONGDFCP/30k1q3bq2TJ09q0aJFGjp0qO688075+fmpfv36mj59um688UaFhoaqQYMGl+2z0qtXL40ZM0Zly5bVgAEDXC7/ZGZmqmrVqrr77rvVuHFjBQQEaNGiRVq9erVeffXVK1pnSRo5cqT+9re/KSEhQYMHD9bbb7+tW265RQ0bNtTDDz+smjVrKjU1VUlJSdq/f7/Wr18vSXryySf1ySefqGPHjnr00Uet2+6rVaumY8eOuZwlGzhwoAYPHqyePXuqY8eOWr9+vRYuXKiwsLBL1la3bl3VqlVLTzzxhA4cOKCgoCDNnj37qvru5HnxxRf17bffqk2bNho0aJDq1aunQ4cOaebMmVq5cqVCQkLcmt+VrmNh1a5dW99++63atm2r2NhYLVmyREFBQVc1z2bNmqlnz56aNGmSjh49at12v23bNklXfpMAcEEeursNKHXyblfOG3x9fU1kZKTp2LGjmTx5ssnIyCjwnovddn+xW9dTU1NNXFyciYqKMj4+PiYyMtK0b9/evPfeey7tTp06Zf7xj3+Y6Ohoq93dd99tdu7cabVZtWqVadasmfH19S30Lczbt2+31m/lypUu07KysszIkSNN48aNTWBgoPH39zeNGzc2U6ZMuex8L/XIgpycHFOrVi1Tq1Ytc+7cOWOMMTt37jR9+vQxkZGRxsfHx1SpUsV07drVzJo1y+W969atM61btzZOp9NUrVrVTJgwwbzxxhtGkklJSXFZxqhRo0xYWJgpV66ciY2NNTt27CjUbfe//vqr6dChgwkICDBhYWHm4YcfNuvXrzeSzLRp06x2ffv2Nf7+/pf9LPLbs2eP6dOnj/X4gZo1a5q4uDiTlZV1yc/tQnUWdh3dmWf+2+7z/PjjjyYwMNDceuut1iMJ/vj9yvveHzlyxOW9ecvO/0iEkydPmri4OBMaGmoCAgJM9+7dzdatW40k89JLLxX2owQuy2EMvdcA2Ed8fLz+/e9/68SJExft1IuSLTk5WU2aNNEnn3xS5I9RgH3RhwhAqfXHn884evSoPv74Y91yyy2EoevEH7ehJE2aNEleXl4ufdiAq0UfIgClVkxMjNq2bat69eopNTVV//nPf5SRkaFnn33W06WhkCZOnKg1a9aoXbt28vb21vz58zV//nwNGjRIUVFRni4PpQiXzACUWk8//bRmzZql/fv3y+FwqGnTpho7duxlb9lHyZGYmKjx48fr119/1YkTJ1StWjX17t1b//jHP4r04aMAgQgAANgefYgAAIDtEYgAAIDtcQG2EHJzc3Xw4EEFBgbyIDAAAK4TxhhlZmaqcuXKl/0dQQJRIRw8eJC7GQAAuE7t27dPVatWvWQbAlEhBAYGSjr/gV7to+gBAMC1kZGRoaioKOs4fikEokLIu0wWFBREIAIA4DpTmO4udKoGAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC2RyACAAC25+3pAgDgeuNweLoC+zLG0xWgtOIMEQAAsD0CEQAAsD0CEQAAsD0CEQAAsD0CEQAAsD0CEQAAsD0CEQAAsD2PBqJx48bJ4XC4DHXr1rWmnzlzRnFxcapQoYICAgLUs2dPpaamusxj79696tKli8qVK6eKFStq5MiROnfunEubZcuWqWnTpnI6napdu7YSEhKuxeoBAIDrhMfPEN100006dOiQNaxcudKaNnz4cM2dO1czZ87U8uXLdfDgQd11113W9JycHHXp0kXZ2dlatWqVPvzwQyUkJGjMmDFWm127dqlLly5q166dkpOTFR8fr4EDB2rhwoXXdD0BAEDJ5TDGc8/9HDdunL788kslJycXmJaenq7w8HB99tlnuvvuuyVJW7ZsUb169ZSUlKRWrVpp/vz56tq1qw4ePKiIiAhJ0rvvvqtRo0bpyJEj8vX11ahRo/T1119r48aN1rzvvfdepaWlacGCBYWqMyMjQ8HBwUpPT1dQUNDVrziA6xpPqvYcnlQNd7hz/Pb4GaLt27ercuXKqlmzph544AHt3btXkrRmzRqdPXtWHTp0sNrWrVtX1apVU1JSkiQpKSlJDRs2tMKQJMXGxiojI0ObNm2y2uSfR16bvHlcSFZWljIyMlwGAABQenk0ELVs2VIJCQlasGCB3nnnHe3atUutW7dWZmamUlJS5Ovrq5CQEJf3REREKCUlRZKUkpLiEobypudNu1SbjIwMnT59+oJ1TZgwQcHBwdYQFRVVFKsLAABKKI/+uGvnzp2tfzdq1EgtW7ZU9erVNWPGDPn5+XmsrtGjR2vEiBHW64yMDEIRAAClmMcvmeUXEhKiG2+8UTt27FBkZKSys7OVlpbm0iY1NVWRkZGSpMjIyAJ3neW9vlyboKCgi4Yup9OpoKAglwEAAJReJSoQnThxQjt37lSlSpXUrFkz+fj4aPHixdb0rVu3au/evYqJiZEkxcTEaMOGDTp8+LDVJjExUUFBQapfv77VJv888trkzQMAAMCjgeiJJ57Q8uXLtXv3bq1atUo9evRQmTJldN999yk4OFgDBgzQiBEjtHTpUq1Zs0b9+/dXTEyMWrVqJUnq1KmT6tevr969e2v9+vVauHChnnnmGcXFxcnpdEqSBg8erN9++01PPvmktmzZoilTpmjGjBkaPny4J1cdAACUIB7tQ7R//37dd999Onr0qMLDw3XLLbfohx9+UHh4uCTp9ddfl5eXl3r27KmsrCzFxsZqypQp1vvLlCmjefPmaciQIYqJiZG/v7/69u2r5557zmoTHR2tr7/+WsOHD9fkyZNVtWpVTZ06VbGxsdd8fQEAQMnk0ecQXS94DhGA/HgOkedwxII7rqvnEAEAAHgagQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANheiQlEL730khwOh+Lj461xZ86cUVxcnCpUqKCAgAD17NlTqampLu/bu3evunTponLlyqlixYoaOXKkzp0759Jm2bJlatq0qZxOp2rXrq2EhIRrsEYAAOB6USIC0erVq/Xvf/9bjRo1chk/fPhwzZ07VzNnztTy5ct18OBB3XXXXdb0nJwcdenSRdnZ2Vq1apU+/PBDJSQkaMyYMVabXbt2qUuXLmrXrp2Sk5MVHx+vgQMHauHChdds/QAAQAlnPCwzM9PccMMNJjEx0bRp08Y8/vjjxhhj0tLSjI+Pj5k5c6bVdvPmzUaSSUpKMsYY88033xgvLy+TkpJitXnnnXdMUFCQycrKMsYY8+STT5qbbrrJZZm9evUysbGxha4xPT3dSDLp6elXupoAShGJwVMD4A53jt8eP0MUFxenLl26qEOHDi7j16xZo7Nnz7qMr1u3rqpVq6akpCRJUlJSkho2bKiIiAirTWxsrDIyMrRp0yarzR/nHRsba83jQrKyspSRkeEyAACA0svbkwv//PPPtXbtWq1evbrAtJSUFPn6+iokJMRlfEREhFJSUqw2+cNQ3vS8aZdqk5GRodOnT8vPz6/AsidMmKDx48df8XoBAIDri8fOEO3bt0+PP/64Pv30U5UtW9ZTZVzQ6NGjlZ6ebg379u3zdEkAAKAYeSwQrVmzRocPH1bTpk3l7e0tb29vLV++XG+88Ya8vb0VERGh7OxspaWlubwvNTVVkZGRkqTIyMgCd53lvb5cm6CgoAueHZIkp9OpoKAglwEAAJReHgtE7du314YNG5ScnGwNzZs31wMPPGD928fHR4sXL7bes3XrVu3du1cxMTGSpJiYGG3YsEGHDx+22iQmJiooKEj169e32uSfR16bvHkAAAB4rA9RYGCgGjRo4DLO399fFSpUsMYPGDBAI0aMUGhoqIKCgvToo48qJiZGrVq1kiR16tRJ9evXV+/evTVx4kSlpKTomWeeUVxcnJxOpyRp8ODBeuutt/Tkk0/qoYce0pIlSzRjxgx9/fXX13aFAQBAieXRTtWX8/rrr8vLy0s9e/ZUVlaWYmNjNWXKFGt6mTJlNG/ePA0ZMkQxMTHy9/dX37599dxzz1ltoqOj9fXXX2v48OGaPHmyqlatqqlTpyo2NtYTqwQAAEoghzHGeLqIki4jI0PBwcFKT0+nPxEAORyersC+OGLBHe4cvz3+HCIAAABPIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADbIxABAADb83an8ebNm/X5559rxYoV2rNnj06dOqXw8HA1adJEsbGx6tmzp5xOZ3HVCgAAUCwcxhhzuUZr167Vk08+qZUrV+ovf/mLWrRoocqVK8vPz0/Hjh3Txo0btWLFCmVkZOjJJ59UfHx8qQpGGRkZCg4OVnp6uoKCgjxdDgAPczg8XYF9Xf6IBfw/d47fhTpD1LNnT40cOVKzZs1SSEjIRdslJSVp8uTJevXVV/X000+7VTQAAICnFOoM0dmzZ+Xj41PombrbvqTjDBGA/DhD5DmcIYI73Dl+F6pT9eXCTVpamlvtAQAAShK37zJ7+eWXNX36dOv1PffcowoVKqhKlSpav359kRYHAABwLbgdiN59911FRUVJkhITE5WYmKj58+erc+fOGjlyZJEXCAAAUNzcuu1eklJSUqxANG/ePN1zzz3q1KmTatSooZYtWxZ5gQAAAMXN7TNE5cuX1759+yRJCxYsUIcOHSRJxhjl5OQUbXUAAADXgNtniO666y7df//9uuGGG3T06FF17txZkrRu3TrVrl27yAsEAAAobm4Hotdff101atTQvn37NHHiRAUEBEiSDh06pKFDhxZ5gQAAAMWtUM8hkqQxY8bozjvvVLNmzYq7phKH5xAByI/nEHkOzyGCO4r8OUSStH//fnXu3FlVq1bVkCFDtGDBAmVnZ191sQAAAJ5W6ED0wQcfKCUlRf/9738VGBioxx9/XGFhYerZs6c++ugjHTt2rDjrBAAAKDaFvmR2IZs3b9bcuXP11Vdfac2aNWrRooXuuOMO3XfffapSpUpR1ulRXDIDkB+XzDyHS2ZwhzvH76sKRPkdPnxYc+fO1Zw5c9S6dWs98cQTRTHbEoFABCA/ApHnEIjgDo8EotKMQAQgPwKR53DEgjuKpVO1JC1dulSvvvqqvv/+e0nSv//9b1WrVk3h4eF6+OGHdfr06SuvGgAAwEMK/Ryi999/X0OGDFF0dLT+8Y9/aOzYsXrhhRfUu3dveXl56ZNPPlGFChX00ksvFWe9AAAARa7Ql8waNGigRx55RI8++qgWLFigbt26aerUqerbt68kaebMmRo9erR27NhRrAV7ApfMAOTHJTPP4ZIZ3FEsl8x+++033XHHHZKk2267TQ6HQy1atLCmt2zZ0vqNs8J655131KhRIwUFBSkoKEgxMTGaP3++Nf3MmTOKi4tThQoVFBAQoJ49eyo1NdVlHnv37lWXLl1Urlw5VaxYUSNHjtS5c+dc2ixbtkxNmzaV0+lU7dq1lZCQ4FadAACgdCt0IDpz5oz8/Pys106nU06n0+X1H4PI5VStWlUvvfSS1qxZo59//ll//etfdeedd2rTpk2SpOHDh2vu3LmaOXOmli9froMHD+quu+6y3p+Tk6MuXbooOztbq1at0ocffqiEhASNGTPGarNr1y516dJF7dq1U3JysuLj4zVw4EAtXLjQrVoBAEApZgrJy8vL7Nixw6Snp5u0tDQTGBho1q9fb9LT0016errZtm2b8fLyKuzsLqp8+fJm6tSpJi0tzfj4+JiZM2da0zZv3mwkmaSkJGOMMd98843x8vIyKSkpVpt33nnHBAUFmaysLGOMMU8++aS56aabXJbRq1cvExsbW+ia0tPTjSSTnp5+NasGoJQ4f+GGwRMD4A53jt+FPkNkjNGNN96o8uXLKzQ0VCdOnFCTJk1Uvnx5lS9fXnXq1LmqYJaTk6PPP/9cJ0+eVExMjNasWaOzZ8+qQ4cOVpu6deuqWrVqSkpKkiQlJSWpYcOGioiIsNrExsYqIyPDOsuUlJTkMo+8NnnzuJCsrCxlZGS4DAAAoPQq9F1mS5cuLZYCNmzYoJiYGJ05c0YBAQH64osvVL9+fSUnJ8vX11chISEu7SMiIpSSkiJJSklJcQlDedPzpl2qTUZGhk6fPu1yGTDPhAkTNH78+KJaRQAAUMIVOhC1adOmWAqoU6eOkpOTlZ6erlmzZqlv375avnx5sSyrsEaPHq0RI0ZYrzMyMhQVFeXBigAAQHEqdCAqLr6+vqpdu7YkqVmzZlq9erUmT56sXr16KTs7W2lpaS5niVJTUxUZGSlJioyM1E8//eQyv7y70PK3+eOdaampqQoKCrrg2SGpYIdxAABQuhW6D1GZMmUKNVyt3NxcZWVlqVmzZvLx8dHixYutaVu3btXevXsVExMjSYqJidGGDRt0+PBhq01iYqKCgoJUv359q03+eeS1yZsHAABAoc8QGWNUvXp19e3bV02aNCmShY8ePVqdO3dWtWrVlJmZqc8++0zLli3TwoULFRwcrAEDBmjEiBEKDQ1VUFCQHn30UcXExKhVq1aSpE6dOql+/frq3bu3Jk6cqJSUFD3zzDOKi4uzzvAMHjxYb731lp588kk99NBDWrJkiWbMmKGvv/66SNYBAACUAoW9dW316tVm8ODBJiQkxDRp0sS8+eab5tixY1dxM5wxDz30kKlevbrx9fU14eHhpn379ubbb7+1pp8+fdoMHTrUlC9f3pQrV8706NHDHDp0yGUeu3fvNp07dzZ+fn4mLCzM/P3vfzdnz551abN06VJz8803G19fX1OzZk0zbdo0t+rktnsA+Xn61nM7D4A73Dl+u/1r92fOnNGsWbM0bdo0/fDDD+rWrZsGDBigjh07Fk9iKwH46Q4A+fHTHZ7j3hELdldsv3YvSWXLltWDDz6oxYsXa+PGjTp8+LBuu+02HTt27IoLBgAA8KQrusts//79SkhIUEJCgk6dOqWRI0dy5gQAAFy3Ch2IsrOz9cUXX+g///mPVqxYoc6dO2vSpEnq3LlzkdxdBgAA4CmFDkSVKlVSYGCg+vbtqylTpqhixYqSpJMnT7q040wRAAC43hS6U7WX1/93N3JcoEehMUYOh0M5OTlFV10JQadqAPnRqdpz6FQNd7hz/Pb4b5kBAAB4msd/ywwAAMDTCnXb/R/7CRV1ewAAAE8qVCCqXbu2XnrpJR06dOiibYwxSkxMVOfOnfXGG28UWYEAAADFrVCXzJYtW6ann35a48aNU+PGjdW8eXNVrlxZZcuW1fHjx/Xrr78qKSlJ3t7eGj16tB555JHirhsAAKDIuPXTHXv37tXMmTO1YsUK7dmzR6dPn1ZYWJiaNGmi2NjYUvtMIu4yA5Afd5l5DneZwR3uHL/d/i0zOyIQAciPQOQ5HLHgjmL9LTMAAIDShkAEAABsj0AEAABsj0AEAABsj0AEAABsr1DPIfrll18KPcNGjRpdcTEAAACeUKhAdPPNN8vhcFi/aH8ppfHX7gEAQOlWqEtmu3bt0m+//aZdu3Zp9uzZio6O1pQpU7Ru3TqtW7dOU6ZMUa1atTR79uzirhcAAKDIFeoMUfXq1a1//+1vf9Mbb7yh22+/3RrXqFEjRUVF6dlnn1X37t2LvEgAAIDi5Han6g0bNig6OrrA+OjoaP36669FUhQAAMC15HYgqlevniZMmKDs7GxrXHZ2tiZMmKB69eoVaXEAAADXQqEumeX37rvvqlu3bqpatap1R9kvv/wih8OhuXPnFnmBAAAAxe2Kftz15MmT+vTTT7VlyxZJ588a3X///fL39y/yAksCftwVQH78uKvn8OOucIc7x2+3zxBJkr+/vwYNGnRFxQEAAJQ0V/Sk6o8//li33HKLKleurD179kiSXn/9dX311VdFWhwAAMC14HYgeueddzRixAh17txZx48ftx7EWL58eU2aNKmo6wMAACh2bgeiN998U++//77+8Y9/yNv7/6+4NW/eXBs2bCjS4gAAAK4FtwPRrl271KRJkwLjnU6nTp48WSRFAQAAXEtuB6Lo6GglJycXGL9gwQKeQwQAAK5Lbt9lNmLECMXFxenMmTMyxuinn37Sf//7X02YMEFTp04tjhoBAACKlduBaODAgfLz89MzzzyjU6dO6f7771flypU1efJk3XvvvcVRIwAAQLG6ogcz5jl16pROnDihihUrFmVNJQ4PZgSQHw9m9BwezAh3uHP8vqLnEJ07d06LFi3Sxx9/LD8/P0nSwYMHdeLEiSuZHQAAgEe5fclsz549uu2227R3715lZWWpY8eOCgwM1Msvv6ysrCy9++67xVEnAABAsXH7DNHjjz+u5s2b6/jx49bZIUnq0aOHFi9eXKTFAQAAXAtunyFasWKFVq1aJV9fX5fxNWrU0IEDB4qsMAAAgGvF7TNEubm51s915Ld//34FBgYWSVEAAADXktuBqFOnTi6/WeZwOHTixAmNHTtWt99+e1HWBgAAcE24fdv9/v37FRsbK2OMtm/frubNm2v79u0KCwvTd999Vypvwee2ewD5cdu953DbPdzhzvH7ip5DdO7cOX3++ef65ZdfdOLECTVt2lQPPPCASyfr0oRABCA/ApHnEIjgDneO3253qpYkb29vPfjgg1dUHAAAQElzRYFo69atevPNN7V582ZJUr169TRs2DDVrVu3SIsDAAC4FtzuVD179mw1aNBAa9asUePGjdW4cWOtXbtWDRs21OzZs4ujRgAAgGLldh+iWrVq6YEHHtBzzz3nMn7s2LH65JNPtHPnziItsCSgDxGA/OhD5Dn0IYI7ivW3zA4dOqQ+ffoUGP/ggw/q0KFD7s4OAADA49wORG3bttWKFSsKjF+5cqVat25dJEUBAABcS253qr7jjjs0atQorVmzRq1atZIk/fDDD5o5c6bGjx+vOXPmuLQFAAAo6dzuQ+TlVbiTSg6H44I/8XE9og8RgPzoQ+Q59CGCO4r1OUS5ublXXBgAAEBJ5HYfIgAAgNKm0IEoKSlJ8+bNcxn30UcfKTo6WhUrVtSgQYOUlZVV5AUCAAAUt0IHoueee06bNm2yXm/YsEEDBgxQhw4d9NRTT2nu3LmaMGFCsRQJAABQnAodiJKTk9W+fXvr9eeff66WLVvq/fff14gRI/TGG29oxowZxVIkAABAcSp0IDp+/LgiIiKs18uXL1fnzp2t13/605+0b9++oq0OAADgGih0IIqIiNCuXbskSdnZ2Vq7dq31HCJJyszMlI+PT9FXCAAAUMwKHYhuv/12PfXUU1qxYoVGjx6tcuXKuTyZ+pdfflGtWrXcWviECRP0pz/9SYGBgapYsaK6d++urVu3urQ5c+aM4uLiVKFCBQUEBKhnz55KTU11abN371516dJF5cqVU8WKFTVy5EidO3fOpc2yZcvUtGlTOZ1O1a5dWwkJCW7VCgAASq9CB6Lnn39e3t7eatOmjd5//329//778vX1taZ/8MEH6tSpk1sLX758ueLi4vTDDz8oMTFRZ8+eVadOnXTy5EmrzfDhwzV37lzNnDlTy5cv18GDB3XXXXdZ03NyctSlSxdlZ2dr1apV+vDDD5WQkKAxY8ZYbXbt2qUuXbqoXbt2Sk5OVnx8vAYOHKiFCxe6VS8AACid3H5SdXp6ugICAlSmTBmX8ceOHVNAQIBLSHLXkSNHVLFiRS1fvly33nqr0tPTFR4ers8++0x33323JGnLli2qV6+ekpKS1KpVK82fP19du3bVwYMHrT5O7777rkaNGqUjR47I19dXo0aN0tdff62NGzday7r33nuVlpamBQsWXLYunlQNID+eVO05PKka7ijWX7sPDg4uEIYkKTQ09KrCkHQ+bOXNS5LWrFmjs2fPqkOHDlabunXrqlq1akpKSpJ0/vlIDRs2dOnwHRsbq4yMDOsxAUlJSS7zyGuTN48/ysrKUkZGhssAAABKrxLzpOrc3FzFx8frL3/5ixo0aCBJSklJka+vr0JCQlzaRkREKCUlxWqTPwzlTc+bdqk2GRkZOn36dIFaJkyYoODgYGuIiooqknUEAAAlU4kJRHFxcdq4caM+//xzT5ei0aNHKz093Rp4nAAAAKWb2z/uWhyGDRumefPm6bvvvlPVqlWt8ZGRkcrOzlZaWprLWaLU1FRFRkZabX766SeX+eXdhZa/zR/vTEtNTVVQUJD8/PwK1ON0OuV0Ootk3QAAQMnn0TNExhgNGzZMX3zxhZYsWaLo6GiX6c2aNZOPj48WL15sjdu6dav27t2rmJgYSVJMTIw2bNigw4cPW20SExMVFBSk+vXrW23yzyOvTd48AACAvbl9l1lRGjp0qD777DN99dVXqlOnjjU+ODjYOnMzZMgQffPNN0pISFBQUJAeffRRSdKqVasknb/t/uabb1blypU1ceJEpaSkqHfv3ho4cKBefPFFSedvu2/QoIHi4uL00EMPacmSJXrsscf09ddfKzY29rJ1cpcZgPy4y8xzuMsM7nDr+G08SNIFh2nTplltTp8+bYYOHWrKly9vypUrZ3r06GEOHTrkMp/du3ebzp07Gz8/PxMWFmb+/ve/m7Nnz7q0Wbp0qbn55puNr6+vqVmzpssyLic9Pd1IMunp6VezugBKifOHZQZPDIA73Dl+e/QM0fWCM0QA8uMMkedwxII7ivU5RAAAAKUNgQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANgegQgAANiet6cLAD8U6Un8UCQAQOIMEQAAAIEIAACAQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGzP29MFAABQYjgcnq7Avozx6OI5QwQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGyPQAQAAGzPo4Hou+++U7du3VS5cmU5HA59+eWXLtONMRozZowqVaokPz8/dejQQdu3b3dpc+zYMT3wwAMKCgpSSEiIBgwYoBMnTri0+eWXX9S6dWuVLVtWUVFRmjhxYnGvGgAAuI54NBCdPHlSjRs31ttvv33B6RMnTtQbb7yhd999Vz/++KP8/f0VGxurM2fOWG0eeOABbdq0SYmJiZo3b56+++47DRo0yJqekZGhTp06qXr16lqzZo3+9a9/ady4cXrvvfeKff0AAMB1wpQQkswXX3xhvc7NzTWRkZHmX//6lzUuLS3NOJ1O89///tcYY8yvv/5qJJnVq1dbbebPn28cDoc5cOCAMcaYKVOmmPLly5usrCyrzahRo0ydOnUKXVt6erqRZNLT06909S7p/OM5GTwxAFfC099bOw9s3FI8FAN3jt8ltg/Rrl27lJKSog4dOljjgoOD1bJlSyUlJUmSkpKSFBISoubNm1ttOnToIC8vL/34449Wm1tvvVW+vr5Wm9jYWG3dulXHjx+/4LKzsrKUkZHhMgAAgNKrxAailJQUSVJERITL+IiICGtaSkqKKlas6DLd29tboaGhLm0uNI/8y/ijCRMmKDg42BqioqKufoUAAECJVWIDkSeNHj1a6enp1rBv3z5PlwQAAIpRiQ1EkZGRkqTU1FSX8ampqda0yMhIHT582GX6uXPndOzYMZc2F5pH/mX8kdPpVFBQkMsAAABKrxIbiKKjoxUZGanFixdb4zIyMvTjjz8qJiZGkhQTE6O0tDStWbPGarNkyRLl5uaqZcuWVpvvvvtOZ8+etdokJiaqTp06Kl++/DVaGwAAUJJ5NBCdOHFCycnJSk5OlnS+I3VycrL27t0rh8Oh+Ph4/fOf/9ScOXO0YcMG9enTR5UrV1b37t0lSfXq1dNtt92mhx9+WD/99JO+//57DRs2TPfee68qV64sSbr//vvl6+urAQMGaNOmTZo+fbomT56sESNGeGitAQBAiVMs97kV0tKlS42kAkPfvn2NMedvvX/22WdNRESEcTqdpn379mbr1q0u8zh69Ki57777TEBAgAkKCjL9+/c3mZmZLm3Wr19vbrnlFuN0Ok2VKlXMSy+95Fad3HZfegfgSnj6e2vngY1biodi4M7x23F+++NSMjIyFBwcrPT09GLpT+RwFPksUUh8+3El2Gc9p9j3WTau5xTDxnXn+F1i+xABAABcKwQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABge96eLgAo1RwOT1dgX8Z4ugIA1xHOEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANsjEAEAANuzVSB6++23VaNGDZUtW1YtW7bUTz/95OmSAABACWCbQDR9+nSNGDFCY8eO1dq1a9W4cWPFxsbq8OHDni4NAAB4mG0C0WuvvaaHH35Y/fv3V/369fXuu++qXLly+uCDDzxdGgAA8DBvTxdwLWRnZ2vNmjUaPXq0Nc7Ly0sdOnRQUlJSgfZZWVnKysqyXqenp0uSMjIyir9YXFNs0lKMjVsqsVlLsWLYuHnHbWPMZdvaIhD9/vvvysnJUUREhMv4iIgIbdmypUD7CRMmaPz48QXGR0VFFVuN8IzgYE9XgGLDxi2V2KylWDFu3MzMTAVfZv62CETuGj16tEaMGGG9zs3N1bFjx1ShQgU5HA4PVlayZGRkKCoqSvv27VNQUJCny0ERYtuWXmzb0ontemHGGGVmZqpy5cqXbWuLQBQWFqYyZcooNTXVZXxqaqoiIyMLtHc6nXI6nS7jQkJCirPE61pQUBA7YCnFti292LalE9u1oMudGcpji07Vvr6+atasmRYvXmyNy83N1eLFixUTE+PBygAAQElgizNEkjRixAj17dtXzZs3V4sWLTRp0iSdPHlS/fv393RpAADAw2wTiHr16qUjR45ozJgxSklJ0c0336wFCxYU6GiNwnM6nRo7dmyBy4u4/rFtSy+2benEdr16DlOYe9EAAABKMVv0IQIAALgUAhEAALA9AhEAALA9AhEAALA9AtF1rG3btoqPj5ck1ahRQ5MmTfJoPUWtNK6Tp+zevVsOh0PJycmeLuWC8n+X7aC077uFlZCQwENvL6Ok77uFtWzZMjkcDqWlpXm6lIsiEJUSq1ev1qBBg4p1Gdf6j9cf18nhcOjLL7+8ZsvHtfO///1Pzz//vPXaTiGhNO67wPWIQFRKhIeHq1y5chedfvbs2WtYzdXJzs6WdPl1Kg3y1vV6UBy15s0zNDRUgYGBRT7/60FJ33c9vfySqjTuu9fTOhUHAtF14uTJk+rTp48CAgJUqVIlvfrqqy7T//g/aofDoXfeeUd33HGH/P399cILL0iSvvrqKzVt2lRly5ZVzZo1NX78eJ07d856X1pamh555BFFRESobNmyatCggebNm6dly5apf//+Sk9Pl8PhkMPh0Lhx4y5Y6xNPPKGuXbtarydNmiSHw6EFCxZY42rXrq2pU6dKkvr166fu3bvrhRdeUOXKlVWnTp0C61SjRg1JUo8ePeRwOKzXhVmnkqRt27YaNmyY4uPjFRYWptjYWG3cuFGdO3dWQECAIiIi1Lt3b/3+++/WezIzM/XAAw/I399flSpV0uuvv17gEtOFzp6FhIQoISHhgnXk5ORowIABio6Olp+fn+rUqaPJkye7tLnYdsnvrbfeUoMGDazXX375pRwOh959911rXIcOHfTMM89IksaNG6ebb75ZU6dOVXR0tMqWLWt9Lnnr07ZtW+3Zs0fDhw+3vmt5Vq5cqdatW8vPz09RUVF67LHHdPLkyYt/4CXA9bTvXs3yX3vtNTVs2FD+/v6KiorS0KFDdeLECZd5JyQkqFq1aipXrpx69Oiho0ePXunHes2Vtn1XOv/de/7559WnTx8FBQVZZyovt599/PHHat68uQIDAxUZGan7779fhw8fdpn3N998oxtvvFF+fn5q166ddu/efZlPuAQwuC4MGTLEVKtWzSxatMj88ssvpmvXriYwMNA8/vjjxhhjqlevbl5//XWrvSRTsWJF88EHH5idO3eaPXv2mO+++84EBQWZhIQEs3PnTvPtt9+aGjVqmHHjxhljjMnJyTGtWrUyN910k/n222/Nzp07zdy5c80333xjsrKyzKRJk0xQUJA5dOiQOXTokMnMzLxgrXPmzDHBwcHm3LlzxhhjunfvbsLCwsyoUaOMMcbs37/fSDLbt283xhjTt29fExAQYHr37m02btxoNm7cWGCdDh8+bCSZadOmmUOHDpnDhw8bY8xl16mkadOmjQkICDAjR440W7ZsMT/88IMJDw83o0ePNps3bzZr1641HTt2NO3atbPeM3DgQFO9enWzaNEis2HDBtOjRw+XbW/M+e39xRdfuCwrODjYTJs2zRhjzK5du4wks27dOmOMMdnZ2WbMmDFm9erV5rfffjOffPKJKVeunJk+fbr1/ottl/x++eUX43A4rO0RHx9vwsLCTK9evazllCtXziQmJhpjjBk7dqzx9/c3t912m1m7dq1Zv3699bnkrc/Ro0dN1apVzXPPPWd914wxZseOHcbf39+8/vrrZtu2beb77783TZo0Mf369buyjXGNXE/77pUu3xhjXn/9dbNkyRKza9cus3jxYlOnTh0zZMgQa/oPP/xgvLy8zMsvv2y2bt1qJk+ebEJCQkxwcHCRft7FpbTtu8ac/+4FBQWZV155xezYscMaLref/ec//zHffPON2blzp0lKSjIxMTGmc+fO1vS9e/cap9NpRowYYbZs2WI++eQTExERYSSZ48ePX9kGuAYIRNeBzMxM4+vra2bMmGGNO3r0qPHz87vkH9X4+HiX+bRv3968+OKLLuM+/vhjU6lSJWOMMQsXLjReXl5m69atF6xj2rRphfrjdfz4cePl5WVWr15tcnNzTWhoqJkwYYJp2bKlMcaYTz75xFSpUsVq37dvXxMREWGysrJc5nOhdfrjH47LrVNJ06ZNG9OkSRPr9fPPP286derk0mbfvn1Gktm6davJyMgwPj4+ZubMmdb0tLQ0U65cuav6o3ohcXFxpmfPntbri22X/HJzc02FChWs+m6++WYzYcIEExkZaYwxZuXKlcbHx8ecPHnSGHM+EPn4+FgBKv/nkn99/rjtjTFmwIABZtCgQS7jVqxYYby8vMzp06cvWqMnXW/77pUu/0JmzpxpKlSoYL2+7777zO233+7SplevXtdVICpN+64x57973bt3dxl3JfvZ6tWrjSQraI8ePdrUr1/fpc2oUaNKfCDiktl1YOfOncrOzlbLli2tcaGhoRc9DZqnefPmLq/Xr1+v5557TgEBAdbw8MMP69ChQzp16pSSk5NVtWpV3XjjjYWu7cUXX3SZ3969exUSEqLGjRtr2bJl2rBhg3x9fTVo0CCtW7dOJ06c0PLly9WmTRuX+TRs2FC+vr6FXm5h16kkatasmfXv9evXa+nSpS71161bV9L57f7bb7/p7NmzatGihfWe4ODgy277wnj77bfVrFkzhYeHKyAgQO+995727t3r0ib/dvn0009d6lyxYoUcDoduvfVWLVu2TGlpafr11181dOhQZWVlacuWLVq+fLn+9Kc/ufSRqV69usLDw92ud/369UpISHCpITY2Vrm5udq1a9fVfRjF5Hrbd690+ZK0aNEitW/fXlWqVFFgYKB69+6to0ePWtM3b97s8jlIUkxMTKHrLQlK076b50Lb+nL72Zo1a9StWzdVq1ZNgYGB1t/zvBqu121tmx93tSN/f3+X1ydOnND48eN11113FWhbtmxZ+fn5ub2MwYMH65577rFeV65cWdL56+3Lli2T0+lUmzZtFBoaqnr16mnlypVavny5/v73v1+y1sK63DqVRPnX9cSJE+rWrZtefvnlAu0qVaqkHTt2FGqeDodD5g8/S3ipzrCff/65nnjiCb366quKiYlRYGCg/vWvf+nHH3+8aK133HGHyx+5KlWqSDq/rd977z2tWLFCTZo0UVBQkBWSLhR+r2ZbP/LII3rssccKTKtWrdoVzbOk8uS+eyXL3717t7p27aohQ4bohRdeUGhoqFauXKkBAwYoOzu71NwcUdr23T+2y1uvS+1nJ0+eVGxsrGJjY/Xpp58qPDxce/fuVWxs7HXfKZtAdB2oVauWfHx89OOPP1p/+I8fP65t27YVONhcStOmTbV161bVrl37gtMbNWqk/fv3a9u2bRf8n6avr69ycnJcxoWGhio0NLRA2zZt2uiDDz6Qt7e3brvtNknnD5z//e9/tW3bNrVt27bQdefx8fEpsPzLrVNJ17RpU82ePVs1atSQt3fB3bFmzZry8fHR6tWrrW2fnp6ubdu26dZbb7XahYeH69ChQ9br7du3X/IM2ffff68///nPGjp0qDVu586dl6w1MDDwgneCtWnTRvHx8Zo5c6a1Xdu2batFixbp+++/LxB+C+NC37WmTZvq119/va629fW4717J8tesWaPc3Fy9+uqr8vI6f+FhxowZLm3q1atX4KD9ww8/FGr5JVFp2Hcvtl6X2s82bNigo0eP6qWXXlJUVJQk6eeff3ZpU69ePc2ZM8dl3PWwrblkdh0ICAjQgAEDNHLkSC1ZskQbN25Uv379rD88hTVmzBh99NFHGj9+vDZt2qTNmzfr888/t+4AatOmjW699Vb17NlTiYmJ2rVrl+bPn2/dHVajRg2dOHFCixcv1u+//37JnfbWW29VZmam5s2b53KQ/PTTT1WpUiW3Tu3nqVGjhhYvXqyUlBQdP368UOtU0sXFxenYsWO67777tHr1au3cuVMLFy5U//79lZOTo8DAQPXt21cjR47U0qVLtWnTJg0YMEBeXl4ud1/99a9/1VtvvaV169bp559/1uDBg+Xj43PR5d5www36+eeftXDhQm3btk3PPvusVq9efUXr0KhRI5UvX16fffaZy7b+8ssvlZWVpb/85S9uz7NGjRr67rvvdODAAeuunVGjRmnVqlUaNmyYkpOTtX37dn311VcaNmzYFdV9LVyP++6VLL927do6e/as3nzzTf3222/6+OOPXe40lKTHHntMCxYs0CuvvKLt27frrbfecrnz9HpTGvbdC7ncflatWjX5+vpa23rOnDkuzxCTzp993L59u0aOHKmtW7fqs88+u+hdcyWKpzsxoXAyMzPNgw8+aMqVK2ciIiLMxIkTXTqiFqYDsjHGLFiwwPz5z382fn5+JigoyLRo0cK899571vSjR4+a/v37mwoVKpiyZcuaBg0amHnz5lnTBw8ebCpUqGAkmbFjx16y5saNG1uda/Pm7XA4zL333uvSrm/fvubOO+8s8P4/rtOcOXNM7dq1jbe3t6levXqh16kk+WPnYWOM2bZtm+nRo4cJCQkxfn5+pm7duiY+Pt7k5uYaY4zJyMgw999/vylXrpyJjIw0r732mmnRooV56qmnrHkcOHDAdOrUyfj7+5sbbrjBfPPNN5fsmHnmzBnTr18/ExwcbEJCQsyQIUPMU089ZRo3bmzN82Lb5ULuvPNO4+3tbXWqzMnJMeXLlzetWrVyaTd27FiXZVzsc0lKSjKNGjUyTqfT5P8z9dNPP5mOHTuagIAA4+/vbxo1amReeOGFQtXoKdfbvnuly3/ttddMpUqVjJ+fn4mNjTUfffRRgU60//nPf0zVqlWNn5+f6datm3nllVeuq07VpW3fvdDNC8Zcfj/77LPPTI0aNYzT6TQxMTFmzpw5BTp+z50719SuXds4nU7TunVr88EHH5T4TtUOY/5w8RJAiXby5ElVqVJFr776qgYMGODpcgAUEvtuyUYfIqCEW7dunbZs2aIWLVooPT1dzz33nCTpzjvv9HBlAC6Ffff6QiACrgOvvPKKtm7dKl9fXzVr1kwrVqxQWFiYp8sCcBnsu9cPLpkBAADb4y4zAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABgewQiAABge/8H2ndIEqfwlMQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the results\n",
    "fig, ax = plt.subplots()\n",
    "# ax.bar(speeds.keys(), speeds.values())\n",
    "ax.bar(\n",
    "    \"direct-write\", speeds[\"direct\"][\"write\"], label=\"Direct write\", color=\"b\"\n",
    ")\n",
    "ax.bar(\n",
    "    \"regular-write\",\n",
    "    speeds[\"regular\"][\"write\"],\n",
    "    label=\"Regular write\",\n",
    "    color=\"r\",\n",
    ")\n",
    "ax.bar(\"direct-read\", speeds[\"direct\"][\"read\"], label=\"Direct read\", color=\"b\")\n",
    "ax.bar(\n",
    "    \"regular-read\", speeds[\"regular\"][\"read\"], label=\"Regular read\", color=\"r\"\n",
    ")\n",
    "ax.set_ylabel(\"Speed (MB/s)\")\n",
    "ax.set_title(\"Direct vs Regular chunking\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d729aa-0fe5-4a2e-8bda-1a79c7b2e6f3",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "As we can see, direct chunking offers immediate speed benefits for both writing and reading.  The reason is that we are bypassing the HDF5 pipeline for dealing with chunked data, and we are directly using the Blosc2 library for that.  This is particularly beneficial for large datasets, where the overhead of the HDF5 pipeline can be very significant.\n",
    "\n",
    "Note that, if you are careful, and you use a cframe for serializing the data with the direct chunking method, the resulting HDF5 file can still be decompressed with any HDF5-enabled application, as long as it has the Blosc2 filter available.  Fortunately, the [hdf5plugin](https://hdf5plugin.readthedocs.io/en/stable/#) comes with the Blosc2 filter, so you can use it in your Python applications without any problem.\n",
    "\n",
    "Finally, as the direct chunking allows for more direct interaction with the Blosc2 library, you can experiment with different blockshapes, compressors, filters and compression levels, to find the best configuration for your specific needs.  In particular, as [Blosc2 supports multidimensional double partitioning](https://www.blosc.org/posts/blosc2-ndim-intro/), you may want to try with different blockshapes to see if you can get better performance for your specific data access pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3970fb6e-88eb-487f-8b95-f71d4faf92f0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "/gnu/store/1w5v338qk5m8khcazwclprs3znqp6f7f-python-3.10.7/bin/python3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
